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ABSTRACT 

We investigate the evolution of angular momentum in Smoothed Particle Hydrodynamic 
(SPH) simulations of galaxy formation, paying particular attention to artificial numerical ef- 
fects. We find that a cold gas disc forming in an ambient hot gas halo receives a strong hydro- 
dynamic torque from the hot gas. By splitting the hydrodynamic force into artificial viscosity 
and pressure gradients, we find that the angular momentum transport is caused not by the 
OO ' artificial viscosity but by the pressure gradients. Using simple test simulations of shear flows, 

\Q | we conclude that the pressure gradient-based viscosity can be divided into two components: 

i one due to the noisiness of SPH and the other to ram pressure. The former is problematic even 

^sO ' with very high resolution because increasing resolution does not reduce the noisiness. On the 

, other hand, the ram pressure effect appears only when a cold gas disc or sheet does not contain 

■ enough particles. In such a case, holes form in the disc or sheet, and then ram pressure from 

intra-hole hot gas, causes significant deceleration. In simulations of galactic disc formation, 
. star formation usually decreases the number of cold gas particles, and hole formation leads to 

the fragmentation of the disc. This fragmentation not only induces further angular momentum 
transport, but also affects star formation in the disc. To circumvent these problem, we modify 
the SPH algorithm, decoupling the cold from the hot gas phases, i.e. inhibiting the hydrody- 
namic interaction between cold and hot particles. This, a crude modelling of a multi-phase 
^ \ fluid in SPH cosmological simulations, leads to the formation of smooth extended cold gas 

discs and to better numerical convergence. The decoupling is applicable in so far as the self- 
gravitating gas disc with negligible external pressure is a good approximation for a cold gas 
disc. 
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1 INTRODUCTION 

Understanding the formation of galactic discs is one of the 
most important, unsolved problems in astrophysics. In the cur- 
rently favoured cold dark matter (CDM) cosmological frame- 
work, in which structure builds up hierarchically I Da vis et all 
1985), discs are assumed to form in the potential wells of viri- 
alised dark matter halos t hrough radiative c ooling ( White & Rees 



Dalcanton. Soergel & Summers 
van den Boschl 1200 lh . In this 



1978]: iFall & Efstathioul Il98 
19971 iMo. Mao & White! Il99 
scenario, baryons are required to retain most of the an- 
gular momentum imparted to them by tidal torques in or- 
der for the resulting centrifugally supported discs to have 
realistic sizes. The conservation of gas angular momen- 
tum is an important assumption made in semi-analytic mod- 
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elling of disc for mation iKauffmann. White & Guiderdonill993t 
I Somerville & Primackl Il999t jCole et alj l200Ct iNagashima et alJ 
boOltlOkamoto & Nagashimal2003l) . " 

However, to date, numerical simulations of galaxy formation 
starting from appropriate CDM initial conditions, and allowing 
just radiative cooling of the gas and no star formation or feed- 
back, find that the infalling gas loses too much angular momen- 
tum. This results in discs forming which are much too small (e.g. 
iNavarro & Whitell994lNavarro. Frenk & White! 1995I) . This prob- 
lem is commonly called "the angular momentum problem". The 
angular momentum losses arise during the hierarchical clustering 
process. At early times in a CDM dominated universe, small dense 
dark matter halos form. Radiative cooling is very efficient in these 
objects and a large fraction of the gas cools into them. As these gas 
rich halos merge to form larger halos their incoming orbital angu- 
lar momentum is drained by dynamical friction and exported to the 
dark matter at the outskirts of the new halos. Much of the original 
angular momentum of the gas is lost through these processes by the 
time it reaches the middle and forms a disc. 
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IWeil. Eke & Efstathioul \ 19981) and lEke. Efstathou & Wrightl 
have shown that if cooling is suppressed until the host ha- 
los are well established, then the numerical simulations yield much 
larger discs. Two ideas have been suggested which might prevent 
the early collapse of small proto-galactic gas clouds. One is that 
cooling may be suppressed by feedback due to energy injected 
by stellar winds and supemovae. Simulations invoking very en- 
ergetic feedback have illustrated the possibility of resolvi ng the 
angular momentum problem in this way. ( Thacker & Couchmad 
1200 it ISommer-Larsen. Gotz & Portinarill2002h . The second idea 
that has been suggested to prevent the formation of small proto- 
galactic clouds is to invoke an alternative form of dark mat- 
ter, "warm dark matter," in which cas e the initial d ensity field 
does not have small scale fluctuations 1 Pagels & Primack 1982). 
ISommer-Larsen & Dolgov| 1200 lh and lGovern^to^t*aLn2002) have 
shown that galaxies formed in this model have larger discs and 
smaller bulges than in simulations with CDM. The angular momen- 
tum problem is potentially a strong clue which can help unravel the 
processes of galaxy formation and the complicated star formation 
and feedback processes involved. For this to be possible, however, 
we must be careful to understand the role of any numerical effects 
which may be important in determining the outcome of galaxy for- 
mation simulations. 

Smoothed particle hydrodynamics (SPH) has been w idely 
used t o study galaxy formation (e . g. iKatz. Hemquist & Weinberg ! 
1992; Evrard. Summers & Davis 1994; Navarro. Frenk & White 
19951 ISteinmetz & Navarrd Il999l : iThacker & Couchmad 12000; 
Steinmetz & Navarro 2002) both because its fully Lagrangian na- 
ture is suited to problems that need a wide dynamic range like 
galaxy formation, and because of its simplicity and robustness 
which make it easy to incorporate into iV-body codes. Despite these 
attractive features, there are problems. First of all, most of the SPH 
implementations utilise an artificial viscosity to capture shocks (but 
also llnutsuka1l2002h . This artificial viscosity can introduce nu- 
merical momentum and angular momentum transport , and spuri- 
ous energy dissipation. Indeed, Som mer-Larsen & Dolgovl 1 200 1 ) 
found that the angular momentum of a simulated galaxy increased 
when they used higher resolution. Another problem arises from 
SPH's intrinsic smoothing properties. Since SPH represents a fluid 
element by smoothing over neighbouring particles, it is not well 
suited for treating large density and velocity gradients. This can 
be a serious problem when a cold gas disc forms through radia- 
tive cooling in an ambient hot gas at the virial temperature. In this 
situation, the cold gas disc is much denser than the hot gas and 
generally rotates faster than the ambient hot medium. In addition, 
because star formation can lead to a decrease in the number of par- 
ticles in the disc (~ 90 % of baryonic matter becomes stars in a 
disc), the effective spatial resolution degrades with time, an effect 
which may play an important role at low redshift. 

In this paper we investigate angular momentum transfer from 
a cold gas disc to the hot halo gas and its effect on the simulation 
outcomes. Although alternative techniques based on Eulerian ap- 
proaches coupled with a grid refinement scheme — Adaptive Mesh 
Refinement (AMR) — have been recently implemented in Cosmol- 
ogy lAbel. Bryan & NormanlEoOOl lTevssieJl2002l and references 
therein), we concentrate here on numerical effects in SPH simula- 
tions of galaxy formation because SPH is by far the most widely 
used method in this area. 

The outline of this paper is as follows. A brief description of 
our simulation code is given in Section 2. In Section 3, we carry 
out cosmological simulations of disc formation in a virialised halo, 
and then demonstrate there is angular momentum transfer from the 



cold gas disc to the ambient hot gas. A forensic study is performed 
in Section 4 using simplified simulations to find the source of the 
problem and the dependence on the numerical resolution. In Sec- 
tion 5 we propose that decoupling of cold gas from ambient hot gas 
can avoid the problems found in earlier sections. The effect of the 
decoupling is shown using simplified simulations and cosmological 
simulations. The results are summarised and discussed in Section 
6. 



2 THE CODE 

We use a modified and extende d version of the p arallel TreeSPH 
code GADGET fcpringel. Yoshida & White! 1200 lh . Unless oth- 
erwise specified, we will use the novel formulation of SPH 
that manifestly conserves energy and entropy when appropriate 
(Springel & Hernauist 2002) throughout this paper. As shown in 
Soringel & Hernauist (2002), this formulation greatly reduces nu- 
merical inaccuracies compared to the more commonly used formu- 
lations. 

GADGET emplo ys an artificial viscosity w hich is the shear- 
reduced version ( Bals arall99llSteinmetzll 996) of the "standard" 
Monaghan & Gingoldj |l9 83l) art ificial viscosity. lLombardi et alj 
1 1999) and Thacker et al. 12000) endorsed this form of the artifi- 
cial viscosity. We set the p arameter a that appears in Eq. 27 of 
ISpring el. Yoshid a & Whitel <200ll) to 0.75. We have checked that 
the choice of the value of this parameter hardly affects our results. 

We also modify the gradient of the smoothing kernel from 
the public version of GADGET according to Thomas & Couchman 
1 1992) to overcome the clumping instabilitv iSchussler^^S chmitt 
Il98ll) GADGET adopts the kernel, W, so-called B 2 -spline 
(Monaghan 1985) 
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for ^ u < I , 
for I < u ^ 1, 
for u > 1 , 



(1) 



where r and h are the particle separation and smoothing length, re- 
spectively, and u = r/h. Note that the smoothing kernel is defined 
over the interval [0, h] and not [0, 2h] as is more common. The ker- 
nel gradient vanishes at u = 0, i.e. the pressure gradient forces 
between two close SPH particles vanishes in the limit of small sep- 
aration. We modify the gradient for it ^ | so that even close pairs 
of SPH particles continue to repel each other, but leave the kernel 
itself unchanged, 



dW_ _ dW_ t _ 1\ 16_ 

du du \ 3 / irh 3 ' 



(2) 



Except in cosmological simulations, where we solve for the 
ionisation, we assume a fixed mean molecular weight, fi = 0.59, 
corresponding to a fully ionised gas of primordial composition. We 
use the adiabatic index, 7 = 5/3 throughout this paper. 



3 COSMOLOGICAL SIMULATIONS OF DISC 
FORMATION 

In order to study disc formation and to concentrate on the investiga- 
tion of numerical angular momentum transfer, we wish to avoid the 
physical angular momentum losses which occur during hierarchical 
clustering as much as possible while retaining a degree of realism. 
To meet these requirements, we choose a halo, from a pre-existing 
TV-body simulation, that is known to have a quiet merger history 
- the redshift of the last major merger is larger than 1. In addition 
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the gas is not allowed to cool radiatively until after z = 1. After 
this both cooling and star formation are allowed. This procedure 
completely suppresses the early formation of proto-galactic clouds 
and leads to quiescent gas accretion f or z < 1, with the resu lt that 
a disc forms with a reasonable size <Weil. E ke & Efstathioi ll99l 
Eke, Efstathou & Wrightll200d) . The details of the simulation are 
described in the following subsections. 



3.1 Initial conditions 

The background cosmology that we assume is a low-density flat 
CDM universe (ACDM). This model is currently the favourite 
amongst hierarchical clustering models. We use the following 
choice of the cosmological parameters: = 0.3, h = i?o/100 
km s" 1 Mpc" 1 = 0.7, A = Ap/ (3ffn) = 0-7, and a s = 0.9. The 
baryon density, fi b is set to 0.04 iNetterfi eld et all2002l) . 

To generate our i nitial conditions, w e use the resimulation 
technique introduced by Frenketa We first perform a dark 

matter only simulation in a 35.325/i _1 Mpc periodic cube. On 
this scale, the density fluctuations are still in the linear regime at 
z — 0. Having completed this simulation, we then select a dark 
halo that has a quiet merger history. The halo's mass is about 
1.3 x 10 12 Mq within the sphere having virial overdensity, 
<5 v ir = 337 at z — 0. To make the new initial conditions, the initial 
density field of the parent simulation is recreated and appropriate 
additional short wavelength perturbations are added to the region 
out of which the halo forms. In this region we also place SPH par- 
ticles in a ratio of 1:1 with dark matter particles. The region ex- 
ternal to this was populated with high mass dark matter particles 
whose function is to reproduce the appropriate tidal fields. The ini- 
tial redshift of the simulation is 50. The masses of the SPH and 
high resolution dark matter particles are ~ 2.6 x lO 6 ft _1 M and 
~ 1.7 x 10 7 h~ 1 M Q , respectively. 

The gravitational softening length for the SPH particles is kept 
fixed in comoving coordinates for z > 3 and after this, it is frozen 
in physical units to a value of 0.5 kpc, in terms of the 'equ ivalent' 
Plummer softening given in Sprine el. Yoshida & White! <200ll) . 
The gravitational force obeys the exact r~ 2 law at r > 2.8e. 
The softening lengths for all particles are defined as e p = e sp h x 
(m p /m sp h) a , where m p is the particle mass of a particle 'p' and 
e sp h and m sp h are the softening and mass of the SPH particles. We 
do not allow the smoothing length to become smaller than a mini- 
mum value of h m i n = 1.4e. 



3.2 Cooling and star formation 

For z < 1 the cooling/heating rate and ionisation state of each 
particle are calculated assuming collisional ionisation equilibrium 
and the presence of an evolving but uniform UV background 
I Haardt & Mada u1ll99rj) by using the fitting formula provided by 
Theun set alj 1 19981) . A primordial composition for the gas is as- 
sumed. Inverse Compton cooling is also considered at z < 1 al- 
though the effect is minor. Since we do not include molecular cool- 
ing, the coolest gas typically has a temperature T co id — 10 4 K. 
We define the gas in an overdense region that has a temperature 
T < 3 x 10 4 K as "cold gas." 

Cold gas particles are eligible to form stars when the following 
criteria are satisfied: (i) the gas particle is in a converging flow (V ■ 
Vi < 0) and (ii) the density of the gas particle is above a threshold 
density (pi > pth). We use pth = 5 x 10 -25 g cm -3 . The value that 



we adopt is higher than the typical va lue used in other cosm olog- 
ical simulations, 2 x 10~ J g cm' " 3 <Katz. Weinberg & Hernquist 
1996). This choice allows us to have sufficient col d gas to observe 
the nu merical effect on the cold phase. Note that B uonomo et alJ 
(2000) found that the criterion on the velocity divergence has no 
sizeable effect on their results, and so criterion (i) may not be 
needed. We ignore the Jeans condition, used by some other authors 
(e.g. iKatz. Wein berg &Hernauisllll99d) . for reasons detailed be- 
low. 

The Jeans condition is usually denoted as 



tli/Cs > tdy 



(3) 



where hi, c s , and tdyn = (4-7rGpi)~2 are the smoothing length, 
the sound speed, and the dynamical time of the gas particle %', 
respectively. Since the hi have no direct physical significance in 
the SPH formalism, and depend for example on the particle mass, 
adopting such a Jeans condition, as given above, would introduce 
an unphysical resolution dependence into the simulations. 

In fact this Jeans condition should be regarded as determin- 
ing the resoluti on limit rather than as a star formation criterion. 
Bate & Burker 3 < 19971) have shown that if the minimum resolv- 
able mass ~ 2A r ng bm sp h, where A ng b is the number of neigh- 
bours used in the SPH calculation, becomes larger than the lo- 
cal Jeans mass, ~ G^^p^^c^, artificial fragmentation may oc- 
cur, and real fragmentation will definitely be suppressed. In our 
adopted SPH implementation, the smoothing length, hi, is defined 
as (47r/3)/ij 3 pi = msphA^gb, and one finds combining these re- 
lations that the above resolution limit is equivalent to 

hi I /3\ 3 l 

— <7T2 - t dyn ~n2t dyn . (4) 

C s \7T/ 

It is clear that the condition represented by Eq.|3|is hardly satisfied 
if the condition represented by Eq. [4] is satisfied. We adopt A^gb 
= 40 in this paper, and this gives, for a gas temperature, T = 10 4 
K, that for densities below p max ~ 1.2 x 10~ 25 g cm -3 the SPH 
treatment of the gas is reliable. Note that this is well below our 
threshold density p t h for star formation. Therefore, our results may 
be affected by artificial effects because of insufficient mass resolu- 
tion. 

When a gas particle is eligible to form stars, the star formation 
rate (SFR) for a particle V is 



dp, 
tit 



Pi 

tdyn 



(5) 



where c* is a dimensionless SFR parameter. This formula corre- 
sponds to the Schmidt law that implies an SFR proportional to 
Pi ' a . The value of c* controls the star formation efficiency. The 
physics of star formation is not understood well enough to predict 
the value of this parameter. It is known for spiral galaxies that the 
star formation time-scale is long compared to the dynamical time- 
scale, so to mimic this, the value of c* must be significantly less 
than unity. We assume a value of c» = 0.05 throughout this paper. 
This is consistent with the relat ively small values of c* used in sim- 
ulations of disc formation (e.g.lKatzl 19921 : IWeil. Eke & EfstathiovJ 



ll99gtlThack er & Co uchmarJl2000l) 7For an SPH particle, of mass 
m gas , which is eligible to form a star, of mass m*, the probability 
of this event occurring during a time-step At is given by 



Psf = 



exp 



c*At 

^■dyii 



(6) 



We use two star formation recipes to study the dependence 
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Figure 1. Distribution of stars and cold gas in the galaxy at z = 0. The 
upper 4 panels and the lower 4 panels are for the conversion run and the 
spawning run respectively. The left panels show the stellar distributions and 
the right panels show the cold gas distributions. The greyscale is coded 
by surface density and three dimensional density for stars and cold gas, 
respectively. The z direction is chosen so that the z-axis becomes parallel 
to the stellar angular momentum of each galaxy. 

of the results on the star formation scheme and also the number 
of particles left in a cold gas disc. In the first recipe, a gas par- 
ticle is completely converted into a stellar particle during At so 
that mg as /m* = 1 in Eq.|6| We call a simulation using this star 
formation scheme a "conversion run". The other recipe allows an 
initial SPH particle to spawn up to three stellar particles with mass 
of m sp h/3, where m sp h is the original mass of the SPH particle, 
so that the values possible for m^/m* are 3, 2 or 1. This scheme 
reduces the rapid decrease in the number of cold gas particles in 
the disc due to star formation, and helps counter the large drop in 
the SPH spatial resolution which otherwise occurs as the cold gas 
is used up. The simulation which employs this scheme is dubbed a 
"spawning run." 

3.3 Results 

In Fig. □ we show the distributions of stars and cold gas in the 
conversion and spawning runs at z — 0. The galaxy in the spawning 



Figure 2. Face on views of the distribution of the cold gas. The left and right 
panels show the galaxy in the conversion run and spawning run respectively. 

run has a smoother stellar distribution than in the conversion run as 
expected. Although the galaxy in the conversion run has a stronger 
bar, we cannot find any significant difference in the distribution of 
the stars when we compare the surface density profiles of these 
galaxies. 

However, the morphologies of the cold gas discs are quite dif- 
ferent. The cold gas disc in the conversion run has a core and ring 
structure: the cold gas disc has large holes and most of the cold 
gas particles are found in dense filaments. The cold gas disc in the 
spawning run is perhaps more realistic with spiral arm-like features, 
though a large fraction of the cold gas lives in the arms. One might 
think that the fragmentation is caused by the Toomre instability. 
However, we find that when we calculate the value of Toomre's 
Q-parameter for the azimuthally averaged surface gas density, that 
Q > 1 is satisfied everywhere and at all redshifts. Because of effi- 
cient star formation, the gas surface density never reaches high val- 
ues. Having ruled out the Toomre instability we need to examine 
the time evolution of the gas distributions to understand the phys- 
ical or numerical mechanisms that cause the break up of the cold 
gas discs. 

The distributions of the cold gas particles are shown in Fig. [2] 
for several redshifts in each simulation. We find that the core-ring 
structure is very common, regardless of the star formation scheme. 
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Figure 4. The hydrodynamic torques acting on the cold gas particles in the 
disc. The results for the spawning run are shown. Azimuthally averaged 
specific torques parallel to the angular momentum of the disc are plotted as 
functions of radius. The solid line represents the total hydrodynamic torque. 
This torque is decomposed into the torque from pressure gradient force (dot- 
ted) and that from the artificial viscosity (dashed). 



The spawning galaxy has a larger-sized cold gas disc at lower red- 
shift than the conversion galaxy. This is consistent with the idea 
that there is angular momentum transfer away from cold gas in a 
way which depends on the number of the cold gas particles. 

Since star formation turns low angular momentum gas par- 
ticles into collisionless stellar particles and the gas particles that 
accrete onto the disc later tend to have larger angular momenta, the 
specific angular momentum of the cold gas disc should monotoni- 
cally increase with time. We plot the evolution of the specific angu- 
lar momentum of the cold gas disc in Fig.|5| Here, we consider the 
material in a sphere of radius 20 h^ 1 kpc, which is centred on the 
galaxy centre at each redshift. The specific angular momentum for 
the conversion run shows surprisingly little evolution, being nearly 
constant. In contrast, the angular momentum of the cold gas disc in 
the spawning run increases monotonically except for the last few 
gigayears. This indicates that more angular momentum is lost from 
the cold gas disc in the conversion run compared to the spawning 
run. 

The lower left panel of Fig. [3] shows the integrated cooling 
rates (the mass in stars and cold gas at each redshift) in the galax- 
ies. The figure shows that the cooling in the conversion run is sup- 
pressed relative to the spawning run. This indicates that a signif- 
icant amount of angular momentum from the cold gas disc in the 




Figure 5. The hydrodynamic torque acting on the cold gas particles in the 
disc in the spawning run. The solid line is identical to the solid line in the 
lower panel of Fig. |4] The dotted line shows the torque when we ignore 
the interaction between the cold gas (T < 3 X 10 4 K) and the hot gas 
(T > 10 5 K) particles. 



conversion run has been transfered to the ambient hot gas, and that 
the cooling rate has decreased because the hot gas has been puffed 
up by this "angular momentum feedback." Since the difference in 
the amount of cold gas in the discs between the two simulations is 
small, the difference in the masses of the cooled baryon (i.e. cold 
gas and stars) is mainly due to the difference in the masses of stellar 
discs. By comparing the evolution of the angular momenta of the 
cold gas discs and the number of cold gas particles in the discs (up- 
per right panel of Fig. one might draw a naive conclusion that at 
least 2000 cold gas particles are needed to suppress the angular mo- 
mentum transfer. However, we have to understand the mechanism 
that causes the angular momentum transfer and the fragmentation 
of the cold gas disc before we can reach definite conclusions. 

To this end, we next investigate the hydrodynamic torques act- 
ing on the cold gas particles. In Fig. |4] we plot the azimuthally 
averaged specific hydrodynamic torques acting on the cold gas par- 
ticles for the spawning simulation for 5 radial bins. The negative 
value indicates that the torque spins down the rotation of the disc. 
As expected from Fig. [3] the absolute magnitude of the torque is 
larger at lower redshift. Surprisingly, the torque is dominated not 
by the artificial viscosity but by the contribution to the force due 
to pressure gradients. The total hydrodynamic torques normalised 
by the angular momenta of the cold gas discs (t z /J z ) at z = 
are -0.84 and -0.90 Gyr -1 for the conversion and spawning runs, 
respectively. This means that the hydrodynamic torque can stop the 
rotation completely in only ~ 1 Gyr. While the number of the cold 
gas particles in the spawning run is more than twice as large as 
that in the conversion run, both discs receive comparable torques at 
z = 0. This implies that the strength of the hydrodynamic torque is 
defined by the distribution (morphology) of the cold gas as well as 
the number of the cold gas particles in the disc. 

To know whether there is significant angular momentum trans- 
fer between the cold gas particles at different radii, we calculate the 
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Figure 3. Evolution of the cold gas disc. The upper left, upper right, lower right, and lower left panels show the specific angular momentum of the cold gas 
disc, the number of cold gas particles in the discs, the mass in the stars and cold gas (integrated cooling rate), and the mass of the cold gas disc, respectively, 
as a function of the age of the universe. The solid and dotted lines indicate the conversion and spawning run, respectively. 



instantaneous hydrodynamic torque acting on the cold gas disc ig- 
noring the interaction between the cold (T < 3 x 10 4 K) and hot 
(T > 10 K) phases. Fig.[5]shows the original torque that is iden- 
tical to the solid line in the lower panel of Fig.[4]and the torque ig- 
noring the hot gas. We find that the hydrodynamic torque becomes 
almost at all radii when we decouple instantaneously the cold and 
hot phases. There is no significant transport of angular momentum 
within the cold gas itself. This confirms results from previous stud- 
ies, which insist that the angular momentum transfer in the galactic 
disc itself cannot be a serious problem over a Hubble time when the 
shear-reduced artificial viscosity is adopted (e.g. Steinmetz 1996). 
Unfortunately, all test simulations have been performed in the ab- 
sence of a surrounding (hot) medium. 

Since our results indicate that the artificial viscosity is not very 



important for the angular momentum transfer, the transfer may be 
caused by fragmentation or may cause the fragmentation. We will 
investigate this point using simplified simulations in the next Sec- 
tion. We note that, the gravitational torque acting on the cold gas 
disc is also significant when the cold gas is fragmented. 

4 SIMPLIFIED SIMULATIONS 
4.1 Shear flows 

Several authors have presented shear flow tests using SPH, and their 
results are promising <LombardietaI]|l999t iThacker et afl fcoOO). 
However, these studies have focused primarily on the variation re- 
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Figure 6. The evolution of the velocity profile in the single temperature shear tests. The results from the SPH runs are shown in the three left-hand columns 
((7Vg as , N nsh ) = (32 3 , 40), (64 3 , 40) and (64 3 , 320) from left to right) and the results from the FD code are shown in the right panels. The top and bottom 
panels corresponds the outputs at t = 293, and 440 Myr for the SPH simulation and t = 286, and 461 Myr for the FD simulation, respectively, and an insert 
in the top left panel shows the initial velocity field. 



suiting from using different SPH implementations. Here, we con- 
sider the impact of varying numerical resolution upon the momen- 
tum transferred across sharp velocity gradients, with or without as- 
sociated sharp density gradients. 

4.7.7 A single temperature test 

We use simulations of a periodic cube of side 10 kpc, contain- 
ing 1O 1O M0 of gas, to investigate the SPH transfer of momentum 
across a discontinuity in the velocity field. The gas is all given a 
temperature of 10 6 K. Particles in the central slab with \z\ < 0.3 
kpc are given a velocity of v x = 50 km s , and the remaining 
gas is set up with v x — —50 km s . At this relative velocity of 
100 km s _1 , it takes ~ 100 Myr to cross the box. The self-gravity 
of the gas is ignored for simplicity, and replaced with an external 
potential of the form 



/ 2-kz \ 


- 1 




V I/box / 





where Lbox is the side length of the simulation box. The choice 
of external potential is not very important because the transfer of 
momentum is not greatly affected by the size of the instabilities 
that it suppresses. To generate relaxed initial conditions, we first 
distribute particles using the rejection method assuming hydro- 
static equilibrium to calculate the density. This system is evolved 
without any shear, while damping the particle temperature and 
velocity until a relaxed state is reached. Then the shear is intro- 



duced and the tests are commenced. Simulations have been run 
with (N gas ,N ngh ) = (32 3 ,40), (64 3 ,40) and (64 3 ,320). A sin- 
gle variant of SPH has been used for the runs in this subsection, 
but we do consider some other popular flavours with different sym- 
metrisations in subsection !4. 1.31 

We have also used the three -dimension al Eulerian fixed-g rid 
hydrodynamic code described bv lOuilis. Ibanez & SaeJ il996) to 
provide a comparison. This finite difference (FD) code employs a 
Riemann solver to compute the numerical viscosity, thus removing 
the need for an artificial viscosity, and has already been used to 
simulate gas stripping from a galaxy by the ram pressur e of the 
intracluster medium (ICM) l Ouilis, Moore & Bower 2000) and the 
evolution of a bubble in the ICM iOuilis. Bower & Baloghll200lh . 
We set up the same initial conditions for the shear test on 151 3 cells. 
This number provides ample resolution and ensures that each cell 
hosts only one phase initially. Random noise is added to the density 
in each cell (Ap/p ^ 0.01), otherwise nothing will happen. 

The evolution of the velocity profile for each simulation is 
shown in Fig.[6| For the SPH runs employing only 40 neighbours, 
the velocity shear decays rapidly owing to momentum transfer 
across the shear boundary. Using 320 neighbours, the transfer of 
momentum is significantly suppressed. In the FD simulation, the 
width of the distribution of v x values at a particular z coordinate 
reflects large scale turbulence that is not apparent in the SPH runs. 
The peak velocities are typically larger and the velocities of the lay- 
ers distant from the contact surfaces remain intact even after 5 box 
crossings of evolution. 
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t (Myr) 

Figure 7. Evolution of the mean x velocity of the cold phase in units where 
1 corresponds to the initial velocity and is reached when all shear has 
gone. Solid, dotted, and dashed lines indicate the shear simulations using 
(7V gas , AT ngb ) = (32 3 ,40), (64 3 ,40) and (64 3 ,320), respectively. At 
early times, the momentum loss is caused largely by artificial viscosity. As 
the relative velocity decreases, pressure gradients provide the more impor- 
tant deceleration. 



By its very nature, SPH is not well-suited to solving problems 
involving large discontinuities. In the following subsection we will 
show more vividly how the SPH and FD methods give qualitatively 
different results regarding turbulence. However, it is instructive to 
understand the difference between the SPH runs considered above. 
As we have seen, the hydrodynamical force can be split into two 
components: pressure gradients and artificial viscosity. The latter 
depends upon the relative velocity of the fluids and the size of the 
interacting volume. This boundary layer should be the same for the 
(32 3 , 40) and (64 3 , 320) simulations, because N ngh /N gaB is iden- 
tical. However, the first and third columns in Fig. [5] provide dras- 
tically different results, so we can infer that the pressure gradients 
must be behind this variation. Fig.Qshows at early times that these 
two simulations do lose momentum at a similar rate. The (64 3 , 40) 
run, and its smaller boundary layer, initially slow down less rapidly. 
As the relative velocity of the gases decreases, the dominant force 
leading to deceleration becomes that coming from pressure gra- 
dients. This depends on JV n gb, such that more neighbours yield 
smaller decelerations. Thus, it seems to be noisiness in the SPH 
smoothing of variables which gives rise to these pressure gradients 
and a significant proportion of the deceleration. The late time evo- 
lution in Fig. \7\ shows the reduced deceleration of the (64 3 ,320) 
run relative to the other two. Note that increasing the number of 
neighbours in SPH calculation significantly slows down the simu- 
lation as well as decreasing the mass resolution. llmaeda & Inutsuki] 
(2002) pointed out that density errors in SPH simulations of shear 
flows can be substantially suppressed by treating particle velocity 
and fluid velocity separately. However, their method is quite slow 
and so far only works with a constant smoothing length (Imaeda, 
private communication). 
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Figure 8. Densities (top), pressures (middle), and temperatures (bottom) 
of the particles in the relaxed initial conditions with A' = 32 3 . The target 
density and pressure calculated from the external potential are given by the 
solid lines. 

4.1.2 A two phase gas 

Now that we have seen how SPH behaves when there is a sharp 
velocity gradient, it is worth investigating the more realistic case 
of a cold slab of gas moving relative to a hotter medium. To this 
end, we have recreated initial conditions with the central slab of 
gas having T co id = 10 J K and the remaining gas left at That = 
10 6 K. The sound crossing time for the cold slab (i.e. 0.6 kpc/c s ) 
is ~ 12 Myr, and ~ 57 % of the mass is in the cold phase. This 
time, when creating the initial conditions, only the temperature of 
the cold phase and velocities were damped, although a maximum 
temperature of Thot is also imposed. Simulations were performed 
using 16 3 , 32 3 , and 64 3 particles. 

Fig.[8|shows the densities, pressures, and temperatures of the 
particles in the relaxed initial conditions using N = 32 3 par- 
ticles. It is apparent that the SPH density and pressure deviate 
from the analytical curves near to the boundary. The prese nce of 
features like these is inevitable with SPH (see lPearce et aljfl999t 
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Figure 9. The evolution of the x—z projections of particles in the shear sim- 
ulations using 16 3 SPH particles (left-hand column), 32 3 particles (middle 
column) and gas density in the FD code (right column). 
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Figure 10. The x—y projections of particles at t = 342 (left) and 489 Myr 
(right) in the shear simulation using 16 3 particles. 



iRitchie & Thomasll200ll) . Hot particles near to the cold dense slab 
overestimate their densities and hence pressures. This causes the 
hot gas to expand away from the dense region, adiabatically cool- 
ing in the process, and creating a gap between the two phases which 
is clearly visible in Fig. [5] This figure also shows that the cold slabs 
are divided into layers, the number of which depends upon the size 
of the SPH smoothing length. 

We now switch on the shear flow as before. The bound- 
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Figure 11. Evolution of the mean x velocity of the cold phase in units 
where 1 corresponds to the initial velocity and is reached when all shear 
has gone. Solid, dotted, and dashed lines indicate the shear simulations 
using 64 3 ,32 3 , and 16 3 particles, respectively. The inflection point at 
t ss 350 Gyr for the lowest resolution ran results from hole formation 
in the cold slab and ram pressure becoming important in enhancing the mo- 
mentum loss. 



ary layer is now unstab le to the Kelvin-Helmholtz instability 
lLandau & Lifshitzll987l) . The presence of a fixed gravitational po- 
tential can in principle stabilise long wavelength modes but for our 
configuration this effect can be largely ignored. The left-hand col- 
umn in Fig.|5|shows the evolution of the x—z projections of parti- 
cles in the N = 16 3 simulation. It is apparent at later times that 
the cold slab is breaking apart. Instabilities are clearly visible at 
t = 98 Myr. At t — 342 Myr an underdensity can be seen at 
(x, y) ~ (4.5, 2). This subsequently grows and the cold phase at 
t — 489 Myr no longer looks like a slab. This is confirmed in 
the face-on projections shown in Fig. 1101 Note that while the mor- 
phology of the phases changes, the membership of each phase is 
constant with cold particles remaining cold and hot ones hot. This 
contrasts starkly with the behaviour of the FD simulation, shown 
in the right-hand column of Fig. [5] where the turbulence mixes the 
phases at the boundary between them. 

In the middle column of Fig. [5] we show the particle distri- 
butions in the N = 32 3 simulation. The instability appears at 
t = 49 Myr, and then vanishes quickly. No significant evolution 
is observed in either this run or that with N = 64 3 , and in nei- 
ther case do the two phases mix at all. Rerunning the N = 32 3 
simulation with a cold gas slab of width 0.4 kpc, rather than 0.6 
kpc, does produce holes. These runs lead us to conclude that the 
holes seen here, and also in the simulations of disc formation, were 
formed due to numerical artifacts. To avoid the formation of holes, 
the cold phase must contain enough particles. This is a more im- 
portant consideration than the number of SPH neighbours being 
used. It is quite hard to maintain enough cold gas disc particles to 
prevent hole formation in galaxy formation simulations, especially 
when star formation is included. 

One consequence of hole formation is shown in Fig. \H] where 
the evolution of the mean velocity of the cold phase in the x direc- 
tion is traced for each SPH simulation. Since there is no mixing 
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Table 1. Hydrodynamic acceleration parallel to the x direction acting on 
the cold phase in each SPH simulation at t = 489 Myr. The first column 
indicates the number of particles in the simulation. The second, third, fourth 
and last columns show the total acceleration, the acceleration from the arti- 
ficial viscosity, the acceleration from the pressure gradients, and the accel- 
eration from the ram pressure respectively. The acceleration is normalised 
by the velocity of the cold phase. The unit is Gyr - 1 . 



(2*.) (Sx.) (2*\ (S3l) 

V«x/ total V%/AV V»i/VP Vfx/R 

16 3 -1.9 -0.21 -1.7 -1.0 

32 3 -0.55 -0.35 -0.21 0.035 

64 3 -0.59 -0.32 -0.27 -0.011 



of the two phases, the evolution of the mean velocity is equivalent 
to the evolution of the momentum of the cold slab. The momen- 
tum of all the gas is well conserved, so any change in the velocity 
of the cold phase should be regarded as the result of a momen- 
tum transfer with the hot phase. There is a resolution dependence 
of the size of the initial deceleration when the slab is perturbed by 
the Kelvin-Helmholtz instability. We have confirmed that this de- 
celeration is caused by the pressure gradient force rather than the 
artificial viscosity. After the initial deceleration, the slabs lose their 
momentum at an almost constant rate, where the acceleration from 
pressure gradients and artificial viscosity are both important. For 
the lowest resolution run, an additional feature in the momentum 
evolution can now be seen at t ~ 350 Myr. This corresponds to the 
epoch when the hole forms in this run and ram pressure from the 
hot intra-hole gas leads to an extra deceleration of the cold mate- 
rial. Such an effect is not evident in the simulations using 32 3 and 
64 3 particles. Their results nearly converge, apart from the initial 
deceleration, although presumably not to a realistic answer given 
that the SPH scheme suppresses any phase mixing. 

In TableQ we show the hydrodynamic acceleration parallel to 
the x direction acting on the cold phase in each run at t = 489 Myr. 
As expected, the cold slab in the N — 64 3 run receives weaker de- 
celeration from the artificial viscosity than in the TV = 32 3 run. 
Although the cold phase in the TV = 16 3 run receives the weakest 
deceleration from the artificial viscosity, we ignore it because the 
morphology of the cold phase in this run is no longer a slab. Inter- 
estingly, the deceleration from pressure gradients is greater in the 
TV = 64 3 run than in the TV = 32 3 run. We also estimate the contri- 
bution of ram pressure to the pressure gradient-based deceleration 
by calculating 

_ 1 dpvl 

ORP — a (o) 

p OX 

using an SPH gather formulation. In the TV — 16 3 simulation, half 
of the deceleration comes from ram pressure, while the contribution 
of the ram pressure is negligible in other simulations. It proves that 
the artificial hole formation in the smallest simulation significantly 
enhances the momentum transfer. 

The slabs in TV = 32 3 and 64 3 receive larger acceleration 
from the artificial viscosity than the pressure gradients. We have 
checked, however, that the slabs lose their momentum more quickly 
only by pressure gradients when we switch off the artificial viscos- 
ity, because the lack of the artificial viscosity increases the numer- 
ical diffusion and the noisiness in the particle distribution. 
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Figure 12. The same as Fig. II II but for the different implementations of 
SPH. In each simulation, 32 3 particles are used. The solid, dotted, and dot- 
dashed lines indicate "entropy: conservation", "energy: asymmetric", and 
"energy: geometric" implementations, respectively. The solid line here is 
the same as the dotted line in Fig. 1111 

4.1.3 Other SPH implementations 

While the numerical difficulties caused by sharp boundaries will 
inevitably impact adversely upon all standard implementations of 
SPH, it is useful to be aware of the variation in momentum trans- 
fer rates within this set of algorithms. We now show the results for 
two other SPH implementations. The first one employs an arith- 
metically symmetrised equation of motion and an asymmetric form 
of the energy equation. This was dubbed 'energy: asymmetric' 
in Soringel & Hernauist (2002), produced the best results among 
the conventional SPH fami ly and has been wide ly used in galaxy 
forma t ion <Evraro1ll988l:lRasio & Shapiroll 199 it Navarro & White 
I l993t ISteinmetz & Mullerl 1 19931 iHaltman & Kallanderl 1 19971 
iThackeretafl EoOO). The second SPH variant we choose is one 
that employs the geometric mean to symmetrise the equation of 
motio n and the energy e quation. This formula was proposed by 
Hgmrjuist_&Katj (l989|), an d was called 'energy: geometric' in 
SDi'ingel & Hernauist 1 2002]) w here it produced the worst re sults 
in their tests. In accordance with lSpringel & Hernauist! J2002F) . we 
refer to our default implementation as 'entropy: conservation.' Note 
that we regenerate initial conditions for each run using the appro- 
priate SPH variant to generate relaxed initial conditions. 

In Fig. [H we plot the evolution of the mean velocity of the 
cold slabs in the simulations adopting the above three flavours of 
SPH. The number of the particles in each simulation is 32 3 . If we 
assume that better algorit hms lose less momentum, then we reach 
the same conclusion as lSpringel & Hernauist! 1200 2) despite using 
a completely different test, namely that 'entropy: conservation' is 
the best and 'energy: geometric' is the worst. 

4.2 Disc formation in a rotating sphere 

We now consider the idealised case of disc formation in a virialised 
rotating sphere similar to that studied bv lNavarro & White! i!99l 
and lThacker et alJ<2000l) . The purpose of this exercise is to obtain 
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a quantitative understanding of how the numerical effects studied in 
the preceding subsections impact upon simulations that are closer 
to what we expect occurs when galaxies form. 

The initial conditions are created by placing particles on a cu- 
bic grid with a spherical edge, and perturbing them radially to give 
a density profile of the form p(r) oc r _1 . Their velocities are cho- 
sen so that the sphere will end up in solid-body rotation around the 
z-axis; the initial angular momentum J corresponds to a value of 
the spin parameter, A = J\E\ 1/2 /(GM 2 ' 5 ) ~ 0.1, where E is 
the total energy of the system. The initial radius of the sphere is 
taken to be 100 kpc, and its mass M 1O 12 M , giving a free-fall 
time from the edge of the system of about 524 Myr. A baryon frac- 
tion of 0. 1 is assumed and equal numbers of dark matter and gas 
particles are used. 

In order to prevent the disc from becoming Toomre unsta- 
ble, we impose a high minimum temperature T m j n = 10 5 K. 
This floor is still well below the virial temperature of the system 
T v i r ~ 3 x 10 6 K, and has the additional benefit of softening the 
Jeans condition given in Eq.|4] 

The simulation is allowed to evolve for 1.25 Gyr without cool- 
ing to let the hot gas reach equilibrium in the halo. Then radiative 
cooling is switched on (using the cooling function computed for 
collisional ionisation equilibrium by Sutherland & Dopit3 1 19931) 
assuming a primordial mix of H and He, and fj, = 0.59) and the 
system is followed until t = 8 Gyr. 
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4.2.1 Cooling only simulations 

We first perform a series of cooling only (no star formation) simu- 
lations using TVgas = 1736, 5544, 15408, 28624, and 44442 SPH 
particles. The same gravitational softening, e = 2 kpc, is adopted 
for gas and dark matter particles for the simulation using 2 x 1736 
particles. Softening lengths for other simulations are chosen as 
e = 2 x (1736/iV g as)^ kpc. The SPH smoothing length is allowed 
to decrease to h m i n = in all runs because this choice results in 
better numerical convergence. Here, the relatively high temperature 
floor (10 5 K) obviates the need to impose a minimum smoothing 
length. 

'Cold gas' is defined as gas that with a temperature lower than 
1.3 x 10 5 K. The top panel of Fig. ll3l shows the evolution of the 
cold gas disc mass. Higher resolution allows higher central gas 
densities, so the better resolved runs have higher cooling rates in 
the first ~ 1 Gyr. After this time, the simulations using more than 
~ 15000 gas particles create almost identical disc masses, and that 
with 5544 SPH particles is only a few per cent lower. 

The evolution of the specific angular momenta of the cold gas 
discs, shown in the lower panel of Fig.^| is more varied. Since the 
cooling is not calculated correctly in the simulation with 1736 gas 
particles, we will ignore this case. Among the remaining simula- 
tions, there is a monotonic increase of angular momentum with in- 
creasing resolution and the evolution converges at 7Vg as = 28624. 
In all cases, the specific angular momenta are decreasing functions 
after t = 3 Gyr for all the simulations that resolve the cooling 
adequately, despite the monotonic increase of specific angular mo- 
mentum with radius when the cooling is switched on. As the total 
angular momentum of the gas is well conserved, this implies that 
there is some transfer within the gas. 

In summary, for this particular test, at least 5000 total gas par- 
ticles are needed before the disc mass is well determined, while 
more than 25000 are required to estimate the disc's angular mo- 
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Figure 13. Evolution of the cold gas discs in the cooling only simulations. 
The solid, dotted, dashed, dot-dashed, and triple-dot-dashed lines denote 
simulations that employ N gaB = 1736, 5544, 15408, 28624, and 44442, re- 
spectively. The upper panel shows the mass of the cold gas disc as a function 
of time (integrated cooling rate). The lower panel shows the specific angular 
momentum of the cold gas disc. 

mentum reliably. Only about 60 per cent of these particles actually 
end up in the disc itself. 



4.2.2 Simulations with star formation 

As seen in Section [3] the specific angular momentum of a cold 
gas disc depends strongly on the number of particles it contains. 
Thus, including star formation with an algorithm that decreases the 
number of cold gas particles should exacerbate the loss of angular 
momentum from the cold gas disc. We now present simulations 
similar to those of the previous section, but also employing the 
'conversion' star formation scheme. A lower threshold density of 
pth = 10~ 25 g cm -3 was adopted, because the cold gas disc is 
more diffuse than in the cosmological simulation as a result of the 
higher minimum temperature allowed for the cold gas (10 5 K). We 
only show results from simulations for which a minimum smooth- 
ing length was not imposed, and omit both the lowest resolution 
simulation, in which cooling was not properly followed, and the 
N SBB = 28624 simulation, which is abased by star formation to the 
low resolution family. A simulation using N ga _ s — 130536 is added 
to study resolution effects. 

In Fig. EH we plot the integrated cooling and star formation 
rates, the mass of the cold gas disc, and the specific angular momen- 
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Figure 14. Evolution of the cold gas discs in the simulation with star formation. The integrated cooling rate, the integrated star formation rate, the remaining 
cold gas mass, and the specific angular momentum of the remaining cold gas are plotted in the upper left, upper right, lower left, and lower right panels, 
respectively. The solid, dotted, dashed, and dot-dashed lines indicate the simulations with 7V gas = 5544, 15408, 44442, and 130536, respectively. 



turn of the cold gas disc. The cooling rates are very similar to those 
in the cooling only simulations, so star formation does not greatly 
affect the cooling rate when a halo contains a sufficient number of 
gas particles (-/V gas ^ 5000). At the end of the simulation there is 
a factor of two difference in the cold gas masses as the resolution 
is varied. However, this is a small fraction of the baryonic mass, 
with only 48, 94, 794, and 2192 cold gas particles remaining in 
the iV gas = 5544, 15408, 44442, and 130536 simulations, respec- 
tively. 

As was seen in the cooling only simulations, the decline of the 
specific angular momentum of the cold gas discs starts from t ~ 3 
Gyr when the rapid accretion of cold gas finishes. These decreases 
are much stronger than those in the cooling only simulations, be- 
cause the cold gas discs now contain fewer particles. Despite this, 



the two highest resolution simulations still manage to produce reas- 
suringly similar results. Compared with the cooling only runs, the 
values of the specific angular momenta of the cold gas discs are sig- 
nificantly higher. This is because the low angular momentum cold 
gas is preferentially creamed off and converted into stars. 

The lower panel in Fig. ll5l shows the final specific angular mo- 
menta of the stars as a function of their formation time. The over- 
all shape of these curves reflects the angular momentum evolution 
in the cold gas, albeit shifted downwards by ~ 30 per cent. The 
declining stellar specific angular momentum at late times implies 
that there would be outside-in disc formation, i.e. older stars have 
larger angular momentum. The outside-in disc formation found by 
Sommer-Larsen, Gotz & Portinari ( 2002) might be explained by 
this process. It is crucial to obtain a numerically robust estimate of 



Table 2. Hydrodynamic torques acting on the cold gas discs parallel to 
the angular momenta of cold gas discs at t = 5.5 Gyr. The second, third, 
and fourth column show the total hydrodynamic torques, the hydrodynamic 
torques from the pressure gradient force, and the hydrodynamic torques 
caused by the artificial viscosity. The torques are normalised by the angular 
momenta. The unit is Gyr - 1 . 











V Joold / hydro 


V ^cold J VP 


V Jcold / AV 


15408 (no SF) 


-0.048 


-0.017 


-0.031 


28624 (no SF) 


-0.035 


-0.013 


-0.021 


44442 (no SF) 


-0.033 


-0.015 


-0.017 


15408 (SF) 


-0.19 


-0.11 


-0.079 


44442 (SF) 


-0.069 


-0.028 


-0.042 


130536 (SF) 


-0.067 


-0.044 


-0.023 
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the angular momentum evolution of the cold gas before the stellar 
population distribution in the disc can be reliably studied. 

The evolution of the total stellar specific angular momentum 
is shown in the upper panel of Fig. 1151 These results are dominated 
by the bulk of the stars which form in the first few Gyr after cool- 
ing is switched on. The small differences in the cold gas angular 
momenta for the two highest resolution runs at these early times 
are sufficient to imprint similar sized differences in the final stellar 
angular momenta. It should be noted that the specific angular mo- 
menta of the stellar discs would depend on the resolution even if 
the cold gas discs had exactly the same specific angular momenta. 
The reason for this is that higher resolution enables higher density 
regions to be resolved at large radii. In addition, a shorter gravita- 
tional softening length gives higher gas density for a given surface 
gas density. Consequently, higher resolution allows higher angular 
momentum gas to form stars. In order to achieve numerical conver- 
gence, we should include self-regulated star formation tuned to give 
the surface density of the star formation ra t e as a function of the 
surfac e gas density iGerritsen & Ickdll997tlSpringel & Hemauistl 
120031 see also Yepes et al. 1997; Hultma n & Pharasyn 1999) as 
the observ ations suggest iKennicuttl 1998). ICarraro. Lia & Chiosi 
( 1998) and Buonomo et alJ fcOOd) have pointed out that inclusion 
of feedback processes, for example, Type la and II supernovae, stel- 
lar winds, and ultraviolet radiation from massive stars, has a signifi- 
cant impact on the evolution of model galaxies. Including feedback 
processes is, however, beyond the scope of this paper. 

Fig. 1161 shows the distribution of the cold gas particles in the 
star formation runs with iV gas = 15408 and 130536. In the iV gas = 
15408 simulation, the cold gas disc has some small holes at t = 3 
Gyr. At t = 6 Gyr the holes have become large and finally the cold 
gas disc shrinks to the centre at t = 8 Gyr. The evolution in the 
iVg as = 130536 simulation shows the same trend, but the disc is 
much smoother. At t = 3 Gyr the disc has beautiful spiral arms, 
and then some small holes appear at t = 6 Gyr. The final disc is 
much more extended than that in the lower resolution simulation, 
but it has acquired unphysically large holes. 

Table.|2|displays the hydrodynamic torques acting on the cold 
gas discs at t — 5.5 Gyr, both in the cooling-only and star forma- 
tion series. We find that the cold gas discs in the simulations with 
star formation received much stronger negative torques compared 
with the cooling only runs with the same initial -/V gas . Part of the 
reason for this is that the torque from the artificial viscosity has a 
strong dependence on the number of cold gas particles. However, 
the contribution to the torque from the pressure gradients is also 
significantly larger in the star formation runs. This is due to ram 
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Figure 15. Evolution of the stellar discs angular momentum (upper panel) 
and formation time-angular momentum distributions in the stellar discs at 
t = 8 Gyr (lower panel). The solid, dotted, dashed, and dot-dashed lines 
indicate the simulations with 7Y gas = 5544, 15408, 44442, and 130536, 
respectively. 

pressure caused by the creation of holes once enough gas particles 
have been converted to stars. In the -/V gas = 130536 simulation, this 
torque becomes stronger than that caused by the artificial viscosity. 
Thus, even in this best-resolved simulation, the cold gas disc is de- 
pleted to an extent where resolution-dependent holes are formed 
and resolution-dependent fractions of the angular momentum are 
lost. 



5 DECOUPLING THE COLD PHASE FROM THE HOT 
PHASE 

In the previous sections we have seen that the angular momentum 
transfer from the cold gas to the hot halo gas is mainly caused by 
numerical problems which are intrinsic to the SPH technique. The 
collapsing rotating sphere test in the previous section reveals that 
the problem becomes worse when star formation is included. The 
reason for this is as follows. Star formation decreases the number 
of particles in the cold gas disc causing the disc to becomes thin- 
ner. We had found in the preceding section that a disc is prone to 
develop holes when the number of particles in it is not sufficient. 
Once these holes appear, ram pressure between the hot and cold gas 
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Figure 16. Distribution of cold gas in the star formation runs, N gSjB = 
15408 (left) and 130536 (right). The grey scale is coded by the three di- 
mensional gas density. 

phases becomes important and a rapid loss of angular momentum 
from the cold gas ensues. 

It is interesting therefore to see what happens when we inhibit 
the angular momentum transfer from cold gas to hot gas in simula- 
tions of galaxy formation. This can be achie ved in a dra stic way by 
decoupling the cold and hot phases. lPearce et allll999l) proposed a 
decoupling technique whereby the contribution to the hot gas parti- 
cles' density from the cold gas particles is explicitly ignored. This 
has the effect of suppressing the overcooling of the hot gas that 
otherwise results. In this scheme hot and cold particles nonethe- 
less interact through mutual hydrodynamical forces. Our approach 
is more radical: hot and cold particles do not interact at all except 
through gravity. The cold and hot phases are treated as indepen- 
dent fluids. Our appro ach is quite simil ar to the multi-phase model 
proposed bv lSemelin & Combesl l20ol in which they represents a 
warm phase by SPH particles and the ISM by sticky particles and 
do not allow the ISM to interact with ambient gas hydrodynami- 
cally. While clearly we miss some real physics as a result of such 
decoupling, the physics that we miss cannot be modelled properly 
using SPH anyway as we have shown. In the following subsections 
we show simulations of disc formation adopting this decoupling 
technique. 

5.1 Disc formation in a rotating sphere 

The technique of decoupling is appropriate in situations in which 
the pressure from the hot gas is unimportant in confining the cold 
gas component and when the pressure at the midplane of the cold 
gas disc is sufficiently greater than the external pressure on the con- 



tact surfaces. In such a situation, the self-gravitating gas disc can be 
a good approximation. We now present simulations including star 
formation, which are similar to those discussed in Section 4.2.2 and 
illustrated in Figs 14-16. 

Fig. ll7l shows the evolution of the cold gas disc. The integrated 
cooling rates are almost the same as in Fig. 1141 The stellar mass is 
also similar to that in the normal simulation, but more mass remains 
as cold gas (about a factor of 2). The evolution of the specific an- 
gular momenta of the cold gas discs are quite different from those 
in Fig. 1141 Apart from the lowest resolution simulation, the simula- 
tions no longer show a decrease in the specific angular momentum 
as a function of time. Even the lowest resolution simulation shows 
a decline which starts much later and the disc ends up with much 
larger specific angular momentum than in the normal SPH simu- 
lation. This proves that most of the angular momentum loss from 
the cold gas disc is due to the hydrodynamic interaction with the 
ambient hot gas, which cannot be dealt with adequately by SPH. A 
resolution dependence is still present, because the accreting gas can 
have large velocity shears before it reaches the temperature below 
which it becomes decoupled from the hot gas. We have also inves- 
tigated the torque that causes the decline of the angular momentum 
of the cold gas disc in the lowest resolution simulation. We find 
that most of the negative torque is coming from the gravitational in- 
teraction with stellar and dark matter particles. We conjecture that 
this is because the density distributions of the cold gas disc and the 
stellar disc are not smooth enough to prevent angular momentum 
transfer due to tidal torques. 

In Fig. 1181 we show the face-on views of the cold gas discs 
in the simulations using -/V gas = 15408 and 130536. Both pro- 
duce smooth extended cold gas discs without any holes. The main 
difference due to resolution is in the three-dimensional density of 
cold gas particles and in the ability to resolve spiral arms. These 
are direct consequences of the differences in the spatial resolution 
of gravity and mass resolution in the SPH. 

The evolution of the stellar discs and their age-angular mo- 
mentum distribution at t = 8 Gyr are presented in Fig. 1191 Ex- 
cept for the lowest resolution simulation, the angular momentum 
evolution shows better convergence than in the normal SPH sim- 
ulations. The reason why the stellar disc in the iVg as = 15408 
simulation has larger specific angular momentum than that in the 
-'Vgas = 44442 simulation is as follows. Although the cold phase 
is decoupled from the hot phase, angular momentum transfer to the 
hot gas during accretion before the gas temperature has reached the 
decoupling threshold (1.3 x 10 K; see Section 4.2.1) is still al- 
lowed. Consequently, the lower resolution simulation produces a 
lower angular momentum cold gas disc. Since the local star forma- 
tion rate is a function of density and we impose a threshold den- 
sity for star formation, the low angular momentum of the cold gas 
disc does not always result in a lower angular momentum for the 
stellar disc. Higher angular momentum gas particles are often not 
dense enough to be eligible for star formation. Thus, the angular 
momentum transfer that brings such gas particles to the inner disc 
where the gas particles can form stars sometimes produces a higher 
angular momentum stellar disc. The increasing resolution is also 
likely to increase the angular momentum of the stellar disc, because 
higher mass resolution allows fragmentation to be resolved in less 
dense environments. As a result, the high angular momentum gas 
particles that cannot form stars in a low resolution simulation are 
allowed to form stars if they reside in the dense structures like spiral 
arms that are resolved with higher resolution (see Fig. ll8> . 

The above complex picture also provides a reasonable expla- 
nation for the fact that, at the highest resolution, the stellar disc in 





the normal SPH simulation (Fig 15) has a slightly higher angular 
momentum than in the simulation with decoupling (Fig 19), de- 
spite the fact that the angular momentum of the cold gas disc in 
the normal SPH simulation is lower than in the simulation with de- 
coupling. The angular momentum transfer from the cold gas makes 
more cold gas particles eligible to form stars and the holes in the 
cold gas disc enhance the density of the cold gas disc at large radii. 
Of course, the pressure from the hot gas also enhances the density 
and thus the star formation rate all over the disc, but, as we have 
seen, this process cannot be modelled appropriately by SPH. 

The age-angular momentum distribution exhibits the outside- 
in disc formation feature again, although this is weaker than in the 
standard SPH simulation. This is not surprising because the star 
forming region shrinks when the surface density of the cold gas 
disc is reduced by star formation. However, as we mentioned above, 



higher resolution allows stars to form in the outskirts of the disc 
even when the surface gas density becomes quite low. Hence the 
outside-in feature becomes weaker with increasing resolution. Even 
though the resolution limits for the iV gas = 44442 and 130536 
simulations (Eq.[4) are far below the threshold density for star for- 
mation (pth = 10~ 25 g cm -3 ), the age-angular momentum distri- 
butions do not converge. Self-regulated star formation may remove 
this resolution dependence as we discussed previously. We will test 
this in future work. 



5.2 Cosmological simulations 

Finally, we present cosmological simulations adopting the decou- 
pling of the cold and hot phases. In this case, we allow warm gas 
(3 x 10 4 < T < 5 x 10 5 K) to interact with both the cold and hot 
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Figure 18. The same as Fig. ll6l but with the decoupling. 



phases so that the simulations do not fail when there are too few 
cold gas particles to perform the SPH calculation. All conditions 
are the same as the simulations presented in Section [5] except for 
the decoupling between the hot and cold phases. 

In Fig. 1201 we show the redshift evolution of the cold gas 
discs in the decoupling simulations. As in Section|5| we adopt two 
star formation schemes, "conversion" and "spawning." The distri- 
butions of the cold gas particles are quite different from those in the 
standard SPH simulations (see Fig. Q and Fig.|2j- In the standard 
SPH simulations, most of the cold gas particles are in the filaments 
and the remaining regions are almost empty. Now the galaxies have 
smooth, extended cold gas discs. These gas discs have spiral arms 
instead of filaments. Since the cold gas particles have a temper- 
ature ~ 10 4 K, that is one order of magnitude smaller than the 
lower limit in the idealised simulations of Section 4, the cold gas 
discs are much thiner. This makes the problems encountered in the 
shear tests more serious. Consequently, the cold gas disc will have 
an unphysical morphology because of numerical effects unless we 
decouple the cold gas from the hot halo gas. Note that the smooth 
density distribution of the cold gas disc significantly decreases an- 
gular momentum transfer due to tidal torques as well. 

Fig. 1211 shows the specific angular momenta of the cold gas 
discs, the integrated cooling rates, and the integrated star forma- 
tion rates in the galaxies as functions of time. The results from 
the standard SPH simulations (Fig. 3) are also plotted for refer- 
ence. Now, with decoupling, the specific angular momenta of the 
cold gas discs become monotonically increasing functions of time. 
There is a small difference between the specific angular momenta 
of the cold gas discs in the two decoupling simulations. This might 
be because the cold gas disc in the spawning case is less affected 
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Figure 19. The same as Fig. ll5l but with decoupling. 



by the angular momentum transfer between the cold gas and the 
warm gas, although we cannot find any significant difference in the 
hydrodynamic and tidal torques in these two simulations. To decide 
whether the high specific angular momentum of the cold gas disc 
in the spawning run is the result of the larger number of cold gas 
particles or a side effect of the multi-mass SPH imposed by this star 
formation prescription, we would have to perform another conver- 
sion simulation with a larger number of the SPH particles leaving 
all other conditions unchanged. In this paper we do not perform 
this test because it would require a much larger number of particles 
in order to satisfy the Jeans condition (Eq.[4j beyond the threshold 
density for star formation. 

The cooling rates and the star formation rates are almost the 
same between the two decoupling simulations, and they are higher 
than in the standard SPH simulations. This proves that the numeri- 
cal angular momentum feedback puffs up the hot halo gas and the 
cooling rate is reduced in the standard SPH simulations. Because 
this effect was not observed in the idealised simulations, we con- 
clude that as the temperature of the cold gas disc becomes lower, 
the numerical angular momentum transfer becomes more problem- 
atic. 

To find out how the angular momentum transfer affects the 
stellar disc, we plot the evolution of the specific angular momenta 
of the stellar discs and the age-angular momentum distributions at 
z — in the decoupling and standard SPH simulations in Fig. 1221 
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Figure 20. Face-on views of the distribution of cold gas. the left panels 
show the galaxy in the simulation with decoupling and the "conversion" 
star formation prescription and the right panels for the simulation with de- 
coupling and the "spawning" star formation prescription. The length is in 
units of h~ 1 kpc. 



Interestingly, except for the standard conversion simulation, all the 
other simulations show similar evolution of the specific angular 
momentum of the stellar discs. On the other hand, except for the 
standard spawning simulation, all the simulations produce simi- 
lar results for the age-angular momentum distributions. These may 
be just a coincidence. Since we impose a relatively high density 
threshold as a star formation criterion, only cold gas particles that 
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Figure 21. Evolution of the galaxies. The top, middle, and bottom panels 
present the specific angular momenta of the cold gas discs, the integrated 
cooling rates, and the integrated star formation rates, respectively. The solid 
and dotted lines indicate the decoupling simulations with the conversion and 
spawning star formation recipes, respectively. The results from the standard 
SPH simulations are also shown by the dashed (conversion) and dot-dashed 
(spawning) lines for reference. 



have small angular momenta can form stars and thus the thresh- 
old density determines the angular momenta of the stellar discs. It 
makes the specific angular momentum evolution of the stellar discs 
almost independent of those of the cold gas discs. 

In the standard SPH simulation with the conversion star for- 
mation prescription, the formation of holes in the cold gas disc re- 
sults in the cold gas being swept out to form a dense ring at large ra- 
dius at high redshifts (see the top left panel of Fig.|2j- Gas particles 
having large angular momenta can form high angular momentum 
stars in this ring. On the other hand, in the standard spawning sim- 
ulation, the gas disc has a large ring at low redshift (see the bottom 
right panel of Fig.|2j- From this ring, the high angular momentum 
population stars of age ~ 2 Gyr are born. Since these holes are 
numerical artifacts as we showed in the shear tests, the large angu- 
lar momentum of the stellar disc and the high angular momentum 
population in the standard SPH simulations are numerical artifacts 
as well. If the gas disc is puffed up by some feedback processes 
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Figure 22. Evolution of the specific angular momenta of stellar discs as 
a function of time and the age-angular momentum distributions of stars at 
z = 0. The solid, dotted, dashed, and dot-dashed lines indicate the de- 
coupling simulations with conversion and spawning star formation and the 
standard SPH simulations with conversion and spawning star formation, re- 
spectively. 



and if we choose a lower threshold density for star formation, the 
simulations with and without decoupling would produce different 
results as in the case of the idealised simulations with the tempera- 
ture floor at T = 10 5 K. The idealised simulations and the cosmo- 
logical simulations presented in this paper suggest that the angular 
momentum transfer from the cold gas disc sometimes results in too 
large and sometimes too small angular momentum of the stellar 
disc depending on the resolution and the modelling of the ISM and 
star formation. It thus must depend on the modelling of feedback 
as well. 



6 SUMMARY AND DISCUSSION 

We have found in a cosmological SPH simulation that a disc of 
cold gas breaks up into filaments and there is a significant transfer 
of angular momentum from the cold gas disc to the ambient hot 
halo gas. The dominant contribution to the hydrodynamical torque 
between the cold and hot phases comes from the pressure gradient 
forces and not from the artificial viscosity. The hot gas is puffed up 
by this angular momentum 'feedback' and this, in turn, can affect 
the cooling rate. 

By using simple shear tests, we find that SPH cannot correctly 
solve problems where there are strong shear flows. When a dense 
cold gas sheet is moving in ambient diffuse hot gas, the gas sheet 
receives strong negative acceleration due to pressure gradients and 



artificial viscosity. By varying the number of neighbouring particles 
involved in the SPH smoothing, we find that the deceleration due 
to pressure gradients is related to the noisiness in SPH variables. 
Moreover, if the cold gas sheet does not contain a sufficient number 
of particles, it undergoes hole formation. Ram pressure from the hot 
gas in these holes leads to further momentum transfer. 

The comparison of simulations using SPH and a finite- 
difference Eulerian code clearly exhibits a shortcoming of the SPH 
method: SPH does not generate the turbulent mixing of the fluid 
components in the shear flow tests that is present with the Eulerian 
code. Instead, with SPH the fluid velocities change to damp out the 
shear flow but without any mixing of fluids. This feature of SPH 
only becomes problematic when there are large velocity gradients. 
Since previous test simulations did not consider such a large shear, 
SPH has produced promising results. However, as we have shown 
in this paper, such a large shear can exist in galaxy formation prob- 
lems where a rotationally supported cold gas disc forms within hot 
halo gas that is mainly supported by thermal pressure. 

The idealised simulations of disc formation in rotating hot 
gas reveal that this problem has a strong dependence on the res- 
olution. Having star formation in a simulation makes the situ- 
ation worse by decreasing the number of cold gas particles. If 
not enough particles are used, the angular momentum of the cold 
gas disc steeply declines with time. Consequently, the stellar disc 
forms from the outside to the inside. This process may explain the 
outside-in form ation of disc galaxies fou nd in the simulations by 
Som mer-Larsen. Gotz & PortinariH2002t) . 

One way to avoid these problems, is to decouple the cold 
and hot gas phases. In the simulations with decoupling, the spe- 
cific angular momenta of the cold gas discs become monotonically 
increasing functions of time as expected on theoretical grounds. 
The main difference from the standard SPH simulations is seen in 
the morphology of cold gas discs. The simulations with decou- 
pling produce smooth extended cold gas discs that do not have 
any holes. These simulations also reveal that the numerical angu- 
lar momentum transfer sometimes increases the specific angular 
momenta of stellar discs and sometimes decreases them, while it 
always decreases the specific angular momenta of cold gas discs. 
This strange feature is caused by the density dependence of star for- 
mation. Hence, the way in which this problem affects stellar discs 
is highly dependent on resolution and on the modelling of subgrid 
physics in the interstellar medium. 

We do not believe that the angular momentum transfer that 
we have identified here is the main explanation for the angular mo- 
mentum problem in galaxy formation simulations, because the dif- 
ference in the specific angular momenta of stellar discs between 
simulations with and without decoupling is not that large, although 
it seems to depend strongly on the modelling of subgrid physics 
such as star formation. At higher resolution, the numerical break- 
ing up of cold gas discs is likely to increase the specific angular 
momentum of the stellar discs by enhancing the cold gas density 
and hence star formation rate in the outskirts of the discs. Anyway, 
as long as the cold gas disc suffers spurious angular momentum 
transfer and has a strange morphology, the properties of the result- 
ing stellar disc are quite unreliable. 

The angular momentum transfer rakes up all the cold gas into 
the star forming region and spurious fragmentation of the disc in- 
duces quite effective star formation in these dense filaments. Con- 
sequently, the problem can be much more serious when one investi- 
gates the details of simulated galaxies like the cold gas distribution, 
the hot gas distribution, the distribution of stellar populations, and 
all observables related to them. The problem may also affect sim- 
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ulations of elliptical galaxies. After the galaxy has exhausted most 
of its cold gas in, for example, a starburst, newly accreting cold gas 
is quickly supplied to the centre regardless of its angular momen- 
tum (newly accreting cold gas implies the existence of halo gas). 
The decoupling of the cold and hot gas phases that we have intro- 
duced offers the opportunity to investigate the detailed structure of 
galaxies by avoiding this spurious angular momentum transfer. 

However, it seems clear that complex physical processes tak- 
ing place in the interstellar medium and in the hot halo gas 
must play a key role in galaxy formation. A code that can solve 
problems involving large shear motions, together with the abil- 
ity to treat a wide dynamic range, is required to study these pro- 
cesses in detail. AMR is an obvious candidate, although it has 
not yet been widely used in this subject, and thus still needs 
substantial testing. On the other hand, substantial refinements 
of SPH have been introduced recen tly which c ould prove use- 
ful in this context, for example Monaghan 1 1989, redefined par- 
ticle velocity by locally ave raged velocity'). lOwen et alj 1 1998, 
tensorial smooth ing kernel s ). iRitchie & Thomas (2001, smoothed 
pressure SPH). Ilnutsukd fcOOl SPH with Riemann solver), 
Imaeda & Inutsuk aH2002l consi s tent pa rticle velocity with fluid ve- 
locity), Kitsionas & Whitworth (2002, adaptive mass resolution), 
and Sprinsel & Hernauist 1 2002, conservation of both entropy and 
energy). We encourage colleagues who have developed and im- 
plemented these refinements to perfor m the shear tests pr esented 
in this paper. The implementation by llmaeda & Inutsuk al (2002) 
is particularly interesting since it substantially suppresses spurious 
density errors in SPH calculations of shear flows. 

The phase decoupling technique that we have introduced may 
be regarded as a crude way of modelling a multi-phase fluid in cos- 
mological simulations and seems to produce much better results 
than standard SPH. This decoupling seems a reasonable approxi- 
mation when the disk is self-gravitating and consists of cold gas and 
stars. As the multi-phase structure of the interstellar medium begins 
to be resolved in a simulation, this approximation must break down 
since the external pressure from halo gas plays a role in confining 
the hot phase of the interstellar medium. Modelling these processes 
remains a challenge which may hold the key for realistic simula- 
tions of the formation of galactic discs. 
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